﻿namespace ImageDupFinder
{
    using System;
    using System.Drawing;
    using System.Drawing.Imaging;
    using ImageDupFinder.PluginCore;

    /// <summary>
    /// Number of bytes used by BRIEF descriptor 
    /// </summary>
    public enum BRIEF_Bytes : uint
    {
        /// <summary>
        /// BRIEF 16, 16*8bit = 128 bits
        /// </summary>
        BRIEF16 = 16,

        /// <summary>
        /// BRIEF 32, 32*8bit = 256 bits
        /// </summary>
        BRIEF32 = 32,

        /// <summary>
        /// BRIEF 64, 64*8bit = 512 bits
        /// </summary>
        BRIEF64 = 64
    }

    /// <summary>
    /// Class used for calculation of "BRIEF" image feature (Binary Robust Independent Elementary Features).
    /// </summary>
    /// <remarks>
    /// BRIEFDescriptor Algorithm:
    /// <para>
    /// 1. Generate predefined number of point pairs. Every point pair value will be stored in memory as single bit.
    ///    To work properly all instances of BRIEFDescriptor need to use the same set of point pairs.
    ///    That's why they are generated in static constructor.
    /// </para><para>
    /// 2. For every image translate point pairs to actual pixels in images.
    ///    This implementation uses one image patch that covers square shaped area located at point (0.1,0.1) 
    ///    and 0.8 width.
    /// </para><para>
    /// 3. To lower impact of noise in image, apply Gaussian blur at every pixel of point pair.
    /// </para><para>
    /// 4. Extract grey values for both pixels in pair.
    /// </para><para>
    /// 5. Calculate BRIEF value of said pair by comparing grey level values.
    ///    Set 1 if first value is higher, otherwise 0.
    /// </para><para>
    /// 6. Store BRIEF value in memory by setting proper bit.
    /// </para><para>
    /// 7. Similarity between two BRIEFDescriptors is number of point pairs with same value.
    ///    This is calculated as Hamming distance of two binary strings. 
    /// </para>
    /// </remarks>
    public class BRIEFDescriptor : IFeature
    {
        /// <summary>
        /// Feature ID.
        /// </summary>
        public const string FeatureName = "BRIEF";

        /// <summary>
        /// Currently used BRIEF format.
        /// </summary>
        private static readonly BRIEF_Bytes currentBRIEF = BRIEF_Bytes.BRIEF64;

        /// <summary>
        /// size of single memory unit
        /// </summary>
        private static readonly uint memoryUnitSize = sizeof(uint);

        /// <summary>
        /// number of memory units
        /// </summary>
        private static readonly uint memoryUnitNum = (uint)currentBRIEF / memoryUnitSize;

        /// <summary>
        /// descriptor length (in bits)
        /// </summary>
        private static readonly uint briefLength = memoryUnitNum * memoryUnitSize * 8;

        /// <summary>
        /// padding for image patch
        /// </summary>
        private static readonly float patchPadding = 0.1f;

        /// <summary>
        /// Smoothing matrix.
        /// </summary>
        private static readonly double[,] kernel = new double[,]
        {
        { 0.00000067, 0.00002292, 0.00019117, 0.00038771, 0.00019117, 0.00002292, 0.00000067 },
        { 0.00002292, 0.00078633, 0.00655965, 0.01330373, 0.00655965, 0.00078633, 0.00002292 },
        { 0.00019117, 0.00655965, 0.05472157, 0.11098164, 0.05472157, 0.00655965, 0.00019117 },
        { 0.00038771, 0.01330373, 0.11098164, 0.22508352, 0.11098164, 0.01330373, 0.00038771 },
        { 0.00019117, 0.00655965, 0.05472157, 0.11098164, 0.05472157, 0.00655965, 0.00019117 },
        { 0.00002292, 0.00078633, 0.00655965, 0.01330373, 0.00655965, 0.00078633, 0.00002292 },
        { 0.00000067, 0.00002292, 0.00019117, 0.00038771, 0.00019117, 0.00002292, 0.00000067 }       
        };

        /// <summary>
        /// point pairs in image patch
        /// </summary>
        private static readonly PointF[] patchA, patchB;

        /// <summary>
        /// Flag. States if feature is properly calculated.
        /// </summary>
        private bool calculated = false;

        /// <summary>
        /// Binary string that stores descriptor bits in (uint) memory units.
        /// </summary>
        private uint[] brief;

        /// <summary>
        /// Initializes static members of the BRIEFDescriptor class.
        /// Static class constructor.
        /// </summary>
        /// <remarks> 
        /// Static constructor generates random point pairs in image patch.
        /// This point pairs are used by all instances of BRIEFDescriptor.
        /// </remarks>
        static BRIEFDescriptor()
        {
            Random random = new Random();
            patchA = new PointF[BRIEFDescriptor.BriefLength];
            patchB = new PointF[BRIEFDescriptor.BriefLength];
            for (int i = 0; i < BRIEFDescriptor.BriefLength; i++)
            {
                patchA[i] = new PointF((float)random.NextDouble(), (float)random.NextDouble());
                patchB[i] = new PointF((float)random.NextDouble(), (float)random.NextDouble());
            }
        }

        /// <summary>
        /// Initializes a new instance of the BRIEFDescriptor class.
        /// </summary>
        /// <param name="image">Image used in feature extraction</param>
        public BRIEFDescriptor(Bitmap image)
        {
            // calculate image padding
            int paddingW = (int)(image.Width * patchPadding);
            int paddingH = (int)(image.Height * patchPadding);

            // calculate smoothing matrix radius
            int radius = kernel.GetLength(0) / 2;

            // prepare array for storing brief bits
            this.brief = new uint[memoryUnitNum];

            BitmapData data = image.LockBits(new Rectangle(0, 0, image.Width, image.Height), ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);
            unsafe
            {
                byte* ptr = (byte*)data.Scan0;
                byte* ptrTmp;
                int kernelLength = kernel.GetLength(0);
                int padding = data.Stride - (kernelLength * 3);
                int byteId = 0;
                int memoryUnitBitnum = (int)memoryUnitSize * 8;
                int startX, startY, stopX, stopY, xr, yr;
                int stride = data.Stride;
                int imageW = image.Width;
                int imageH = image.Height;
                double div, k;

                for (int i = 0; i < briefLength; i++)
                {
                    // translate point pairs to pixels
                    int x1 = (int)(paddingW + (patchA[i].X * (imageW - (2 * paddingW)))),
                        y1 = (int)(paddingH + (patchA[i].Y * (imageH - (2 * paddingH)))),
                        x2 = (int)(paddingW + (patchB[i].X * (imageW - (2 * paddingW)))),
                        y2 = (int)(paddingH + (patchB[i].Y * (imageH - (2 * paddingH))));

                    // ---------------------------
                    // blur  procedure for pixel 1

                    // set pointer at image start
                    ptr = (byte*)data.Scan0;

                    startX = x1 - radius;
                    stopX = startX + kernelLength;
                    startY = y1 - radius;
                    stopY = startY + kernelLength;

                    // move pointer to proper pixel.
                    ptr += (startY * stride) + (startX * 3);

                    double blur1 = 0;
                    div = 0;
                    for (int ky = 0; ky < kernelLength; ky++)
                    {
                        // y offset
                        yr = ky - radius;

                        // if smoothing kernel is outside image skip it or move to next row
                        if (startY + yr < 0)
                        {
                            continue;
                        }

                        if (startY + yr >= imageH)
                        {
                            break;
                        }

                        for (int kx = 0; kx < kernelLength; kx++)
                        {
                            // x offset
                            xr = kx - radius;

                            // if smoothing kernel is outside image skip it
                            if (startX + xr < 0)
                            {
                                continue;
                            }

                            if (stopX + xr < imageW)
                            {
                                ptrTmp = ptr;

                                // move pointer by xr and yr offsets
                                ptrTmp += (yr * stride) + (xr * 3);

                                // get kernel
                                k = kernel[kx, ky];

                                // add to blur value
                                blur1 += ((ptrTmp[0] * .3) + (ptrTmp[1] * .59) + (ptrTmp[2] * .11)) * k;

                                // add kernel to divider
                                div += k;
                            }
                        }
                    }

                    // if divider was less than 1 means some pixels were skipped, update blur value
                    if (div < 1.0 && div != 0)
                    {
                        blur1 /= div;
                    }

                    // ---------------------------
                    // blur  procedure for pixel 2

                    // set pointer at image start
                    ptr = (byte*)data.Scan0;

                    startX = x2 - radius;
                    stopX = startX + kernelLength;
                    startY = y2 - radius;
                    stopY = startY + kernelLength;

                    // move pointer to proper pixel.
                    ptr += (startY * stride) + (startX * 3);

                    double blur2 = 0;
                    div = 0;
                    for (int ky = 0; ky < kernelLength; ky++)
                    {
                        // y offset
                        yr = ky - radius;

                        // if smoothing kernel is outside image skip it or move to next row
                        if (startY + yr < 0)
                        {
                            continue;
                        }

                        if (startY + yr >= imageH)
                        {
                            break;
                        }

                        for (int kx = 0; kx < kernelLength; kx++)
                        {
                            // x offset
                            xr = kx - radius;

                            // if smoothing kernel is outside image skip it
                            if (startX + xr < 0)
                            {
                                continue;
                            }

                            if (stopX + xr < imageW)
                            {
                                ptrTmp = ptr;

                                // move pointer by xr and yr offsets
                                ptrTmp += (yr * stride) + (xr * 3);

                                // get kernel
                                k = kernel[kx, ky];

                                // add to blur value
                                blur2 += ((ptrTmp[0] * .3) + (ptrTmp[1] * .59) + (ptrTmp[2] * .11)) * k;

                                // add kernel to divider
                                div += k;
                            }
                        }
                    }

                    // if divider was less than 1 means some pixels were skipped, update blur value
                    if (div < 1.0 && div != 0)
                    {
                        blur2 /= div;
                    }

                    // ---------------
                    // calcualte BRIEF

                    // get current memory unit
                    byteId = i / memoryUnitBitnum;

                    // move left by 1 to create new position for BRIEF value with its bit set to 0
                    this.brief[byteId] <<= 1;

                    // round blur values
                    double b1 = Math.Round(blur1);
                    double b2 = Math.Round(blur2);

                    // set brief value
                    if (b1 < b2)
                    {
                        this.brief[byteId]++;
                    }
                }
            }

            image.UnlockBits(data);

            this.calculated = true;
        }

        /// <summary>
        /// Gets BRIEF desriptor length in bits.
        /// </summary>
        public static uint BriefLength
        {
            get { return BRIEFDescriptor.briefLength; }
        }

        #region IFeature Members

        /// <summary>
        /// Gets feature name.
        /// </summary>
        public string Name
        {
            get { return FeatureName; }
        }

        /// <summary>
        /// Gets a value indicating whether feature is calculated.
        /// </summary>
        public bool IsCalculated
        {
            get { return this.calculated; }
        }

        /// <summary>
        /// Calculate similarity between two BRIEF features.
        /// </summary>
        /// <param name="other">Second feature used in calcualtion</param>
        /// <returns>Similarity in range 0 - 1, where 0 means "identical" and 1 means "completly different"</returns>
        public float GetSimilarity(IFeature other)
        {
            if (!(other is BRIEFDescriptor) || !other.IsCalculated || !this.IsCalculated)
            {
                return 1;
            }

            BRIEFDescriptor fr = other as BRIEFDescriptor;

            // calculate Hamming distance
            int dist = 0;
            for (int i = 0; i < memoryUnitNum; i++)
            {
                // XOR bytes
                uint val = this.brief[i] ^ fr.brief[i];

                // count bits
                while (val > 0)
                {
                    ++dist;
                    val &= val - 1;
                }
            }

            // set lowest posible similarity to 25% of matching point pairs
            return ((float)dist / briefLength) * 1.5f;
        }

        #endregion
    }
}
